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Abstract 

The linear response of two-dimensional amorphous elastic bodies to an external delta force is 
determined in analogy with recent experiments on granular aggregates. For the generated forces, 
stress and displacement fields, we find strong relative fluctuations of order one close to the source, 
which, however, average out readily to the classical predictions of isotropic continuum elasticity. 
The stress fluctuations decay (essentially) exponentially with distance from the source. Only be- 
yond a surprisingly large distance, 6 « 30 interatomic distances, self- averaging dominates, and the 
quenched disorder becomes irrelevant for the response of an individual conflguration. We argue that 
this self-averaging length b sets also the lower wavelength bound for the applicability of classical 
eigenfrequency calculations. Particular attention is paid to the displacements of the source, allow- 
ing a direct measurement of the local rigidity. The algebraic correlations of these displacements 
demonstrate the existence of domains of slightly different rigidity without, however, revealing a 
characteristic length scale, at least not for the system sizes we are able to probe. 
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I. INTRODUCTION. 



The recent years have seen a tremendous effort to determine the response of granular 



matter subject topoint (delta) sources as indicated in Fig.lU^a). These theoretical 

experimental ja Q, ^,0,11^1 and computational [ill studies have been motivated by the 
desire to understand the static properties of, say, a humble sandpile — to quote an important 
paradigmatic example It has been argued that these aggregates formed under gravity 
as external driving force — alongside with other special "solids" such as jammed colloids, 
emulsions or foams — may not necessarily be described as classical elastic or elastoplastic 
continuum bodies 0,0,0] • Hence, the interest to determine experimentally and by computer 
simulation the linear and quasistatic response to a localized incremental force, in order to 
distinguish between the different models proposed. In a nutshell, stress distributions below 
the source rather close to classical elasticity predictions have been found for standard sand, 
although minor differences seem to appear in the distribution tails 0]. This has prompted 
the more recent focus on the fact that these systems are typically composed of a small 

n n 

number of constituents p], and on the paramount role of the quenched disorder pj. 

In this paper, the point source response problem is carried over to a definitely much 
simpler disordered model system, the two-dimensional amorphous solid formed by quenching 
a Lennard- Jones fluid. It is well known for amorphous materials such as metallic, organic 
or mineral glasses, that their mechanical prope rties are quite different from those of the 
corresponding crystals at the same density They are characterized by a large decrease 

in both the apparent shear and Young's moduli, and a large increase of the yield stress 



associated with a localization of the plastic deformation ll2L These properties have 

nnnrT 

been interpreted in terms of local rearrangements [U, [151, [161, [ITJ due to the heterogeneity 
of the microscopic structure. But these rearrangements have never been identified clearly. 
Particularl y, l ike in granular materials, the role of the quenched stresses is actually a matter 
of debate jlS, Q, as well as the role of local heterogeneities in the elastic constants of 
the materials. One way to answer those questions experimentally is to perform nanoscale 
indentation that is to study the response to a point force. 

As schematically illustrated in Fig.[T](b) we study stress and displacement fields generated 
by an external force acting on the Lennard- Jones beads contained within a small disk. 
Snapshots of the incremental stresses and displacement fields presented in Figs. 121 01 and IH 
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show a rather noisy response. Since these systems behave clearly as classical elastic bodies 
— provided sufficiently large wave lengths and small forces are probed — they provide 
important reference systems, the results from granular matter may be compared with. The 
linear response to point source is equivalent to the (noisy and position dependent) Green's 
function. To our knowledge this study presents the ffist systematic computational study of 
this function, extending some aspects discussed only recently jsf. Specifically, we investigate 
how this noise affects the averages compared to classical continuum theory, and how the 
distributions get narrower with increasing distance from the source, due to self-averaging. 
The spatial correlations of the responses of close sources are studied in order to verify whether 
domains of different rigidity exist as has been argued recently liol . l2fll | . 

This work is in fact the natural sequel of our study 22, where the applicability 
of classical elastic continuum theory on small length scales has been tested by comparing 
with theory the low end of the eigenfrequency spectrum obtained by diagonalization of 
the dynamical matrix. We found that only for system sizes and wave lengths larger than 
a characteristic wave length, ,^ ^ 30 interatomic distances, the eigenfrequencies show the 
degeneracy predicted for a classical isotropic and homogeneous body. This surprising large 
lower limit for classical continuum theory is also seen in the characterization of the non-affine 
field generated by macroscopic deformations (shear or elongation) . Only after coarse-graining 
over distances of order ^ does the non-affine response become negligible. This does in turn 
explain why modes associated with smaller wave lengths do not behave as predicted from 
an approach formulated in terms of affine displacement fields. 

In this paper, we describe ffist briefly some technical points related to the initial samples, 
the computational methods and measurements. In the subsequent Sections ITTTl IIVI and fVl we 
present our numerical results for stress and displacement fields, and their distributions. We 
have regrouped our results following the three different boundary conditions investigated. 
In the ffist section, we demonstrate that the self-averaging is characterized by a length scale 
similar to the critical wave length ^ from our previous study. In the latter two sections 
we analyze the source displacements and their correlations. Our results are summarized 
in Sec. IVII The analytical predictions from classical elasticity theory are outlined in the 
appendix. 
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II. COMPUTATIONAL TECHNICALITIES. 



The initial configurations and their preparation have been to some extent described in 
[3]. Of relevance here is a large ensemble of 16 independent configurations containing each 
10,000 Lennard- Jones particles quenched from T = 1 down to zero temperature following 



a fixed protoco 
dient methods 



using standard molecular dynamics, steepest descent and conjugate gra- 



2J|. Note that, while the particle mass m is strictly monodisperse, we use 



sufficiently polydisperse particle diameters (uniformly distributed between 0.8 and 1.2) to 
prevent crystalline order. The linear size of the periodic boxes is L = 104, the corresponding 
volume fraction 0.925. The mean pressure P = 0.25 was chosen to be close to zero. The two 
Lame coefficients [3,0, l^]; ^ ~ 39.5 and /i ~ 11.7, have been measured directly using 
Hooke's law by applying macroscopic elongation and (pure) shear to the simulation box. We 
recall that the associated Poisson ratio u = X/{X + 2fi) 2/3 is larger than 1/2 which is 
permissible in a 2D solid. Here as everywhere below we have naturally given the numerical 
values in Lennard- Jones units. 

It has been carefully checked that the initial configurations and their monomers are indeed 
at mechanical equilibrium, i.e. are sitting in (local) minima of the energy landscape. The 
linear response to a small external force or imposed displacement can, hence, be described 
by means of the (2A^) x (2A^) dynamical matrix M whose elements depend on the frozen 
tensions ( "quenched stresses" ) and stiffnesses of the links between interacting beads 

H- In 

principle, it is straightforward to solve numerically the linear equations M-U_ = K- Here, 
F_ and U_ are the 2A^-dimensional force and displacement fields respectively containing the 
imposed external body forces and displacements. Since we are considering very large systems 
and standard linear equation solver being of order A^^ we have mainly used fSec. IIIII and lIV|) 
direct steepest descent and/or conjugate gradient methods which are in this case (where the 
neighbor contact lists remains constant) of order A^. The advantage of the direct methods 
is also that they allow to probe the non-linear response regime. We have checked that both 
methods yield the same results for sufficiently small external forces. For comparison, we 
present in Sec. E] results obtained directly from the dynamical matrix. 

In all shown in Fig. ^b), we apply a localized external force of / /rio to all 

the Uq beads contained in small source disks of fixed diameter D. The center of the disk 
refers naturally to the origin of our coordinate system {x,y). The special limit with sources 
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containing only one bead (no = I, D ^ 1) will be used in Sec. El Obviously, the response 
becomes locally less noisy with increasing source size. As we are interested in disorder on 
distances larger than the typical particle distance we have also distributed the source over 
more than one bead. Most of the results reported in Sections IIIII and IIVI are for D = 4 
corresponding to (no) ~ 12 beads. All the source forces considered here point vertically 
downwards. It turns out that an applied force of order one per bead is sufficiently small to 
ensure linear elastic response for the direct minimization methods. (See Fig. El below.) The 
averages are taken over different disk positions in the same configuration, and also over the 
configuration ensemble. 

For mechanical stability, we have either imposed a compensation force of —f^/N on all 
beads or fixed the positions of certain beads, as shown in Fig. ^b). The first method has the 
advantage of being free of any fixed boundary layer making it possible to use the full initial 
periodic box. Care has to be taken however, in this case, for numerical reasons, because 
small drifts of the system cannot be completely avoided. The displacement fields must thus 
be considered in the center of mass frame. Sec. IIVI presents results averaged over nearly 4000 
linear responses obtained with this boundary condition. 

Most of the work presented in this paper f Sec. IIIII lV|) uses instead fixed beads to compen- 
sate the source force. Either we fix all beads in a horizontal layer with \y\ > h and h < L/2 
(Sec. IIII|) . or all beads which are beyond a given number of topological layers around the 
source particles (Sec. IV]) . 

It is well known for elastic bodies in two dimensions that stresses and strains far from 
both source and boundary decrease inversely with the source distance r, i.e. the displacement 
field varies logarithmically. Obviously, the response depends generally on the imposed bound- 
ary conditions. Details of the rigorous analytical treatment, exemplified for the boundary 
conditions studied in the next section, are outlined in the appendix. 

III. RESULTS FOR FIXED TOP AND BOTTOM LAYERS. 

The two snapshots of the forces and displacement fields depicted in Fig. |21 and El show 
the response fields obtained in a small system of linear size L = 32.8 containing only = 
1000 beads, but at the same volume fraction and pressure as the larger samples studied 
quantitatively below. The strength of the forces between beads are represented in Fig. El 
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by the width of the hnes repulsive (tensile) forces being black (gray). Only the incremental 
forces 6f due to the source are given, i.e. the rather strong quenched or residual forces of the 
amorphous body have been subtracted. The force chains visible resemble strongly the ones 
known from granular matter jfjO,^^, although our system is certainly a classical isotropic 
elastic body at large distances 2dl\ . 

The displacement field 6u = u — Ucet indicated by the arrows in Fig. El has been obtained 
by substracting from the total displacement field the displacement field Ucet calculated for 
standard continuum elasticity theory (GET) j25,|2^ following the prescription indicated in 
the appendix. In other words, 6u depicts the noisy response due to the quenched disorder. 
In order to do a comparison of both snapshots, we show in Fig. |3] the residual displacement 
field 6u, and the noise component of the local incremental stress on each particle. In order 
to obtain the noise component, we have substracted the stress calculated with standard 
continuum elasticity theory, and the quenched stresses, from the total stress in the presence 
of a source. The total stress has been calculated here on each particle, using the standard 
Kirkwood definition ^, The noisy part of the incremental stress is then represented 
by an ellipse centered on the particle, whose large principal axis is proportional to the 
largest eigenvalue, and whose small principal axis is proportional to the smallest eigenvalue 
of the residual stress tensor. The directions of these axes give thus the main directions of 
the incremental stresses. The snapshot of Fig. |3] shows clearly that 6u is corrrelated to the 
local incremental stress. To get more quantitative results, we have drawn in the inset of 
Fig- El the distribution of the angles 6 between the residual displacements 6u, and the main 
direction of the incremental stresses (the direction associated to the largest eigenvalue). 
We show a peak for zero angle, with a broad distribution (linear with 6). The residual 
displacement field thus reveals a clear tendency to align with the main direction of the 
incremental stresses. On larger distances, however, we see a vortex like structure for 6u 
similar to the structure revealed by the non-affine displacement field under macroscopic 
strain found in |2S]- The reason for this can be easily understood for the latter case where 
the pressure must become macroscopically constant, and with it the particle density as well. 
This generates the "backfiow" of the non-affine displacement, just like in a uncompressible 
fluid. We recall that the continuum displacement field Ucet for the point source problem fiows 
also back, but on a distance L/2 given by the system size. The size of the vortices measured 
in 6u, however, does not depend on the system size. Note that, from one configuration to the 
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next, the vortices are not located at the same place, and that they disappear after averaging 
the displacement field over many configurations. Moreover, these vortices are not due to 
the natural discretization of our system: they would disappear if the spatial distribution of 
atoms were ordered, as can be infered from the direct computation of the Green function 
(see for example Ref. H)- Note finally that the 3D case is now under study. We would not 
be surprised if the size of the structure involved in the local rearrangements of the atoms 
were smaller in the 3D case, due to the minor effect of disorder and the smaller range of 
elasticity 28|. However, systems with a very large number of atoms have to be studied in 
this case j22| to fit with the continuum limit. 

Fig. shows the vertical normal stress (Tyy generated by one source of diameter D, at 
a distance ?/ = 50 below the source, i.e. just above the fixed beads of the bottom layer. 
As in the snapshot Fig. |21 only the incremental stresses due to the source are shown here. 
To make comparison between the sources of different strengths, the total vertical stress has 
been normalized to 1. 

The stress tensor has been measured, as everywhere in the following, by means of the 
virial definition Q averaged over all beads contained in small rectangular volume elements 
of width 5 and height 3 centered at (x, y). Adopting in this work the sign convention usual 
in granular matter, compressive stresses are taken as positive, i.e. have the same sign as the 
pressure. The size and the aspect ratio of the volume elements were chosen for convenience. 
A typical volume element contains 14 beads, and averages over about 100 interactions which 
takes out some of the noise. On the other hand, it remains small enough to achieve a 
good spatial resolution. Note that a given interaction may contribute to two neighboring 
volume elements. Data points corresponding to two such elements are therefore statistically 
correlated, and the curves appear slightly smoother as they would otherwise. 

Two additional points have to be made here. First, the responses compare already quite 
well with the analytical prediction (bold line) albeit they are not averaged over different 
realizations and despite the fact that the (not given) snapshot of the forces still looks quite 
noisy. This is obviously to be expected for large distances from the source as the response in 
an elastic body should self-average over the noise. While we shall make this more quantitative 
in a moment. Fig. El shows clearly that a distance of order ?/ ^ 50 yields a reasonable — 
although not perfect — self-averaged response. This confirms our finding in that systems 
of size L = 104 show accurately the lowest eigenmodes and can therefore be regarded as free 
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of finite size effects. This motivated our choice of this system size. Note that the continuous 
response and the response averaged over many configurations (bold hne in the Fig.|3J coincide 
at this distance from the point source (see also Fig. IHlon this point). 

Second, we note in Fig. El that the responses are identical for all systems where the forces 
per bead remain of order one or lower. This appears to be independent of the disk diameters 
D despite the additional beads charged for larger disks. Apparently, these differences at the 
source are screened at ?/ = 50 ^ Z^. The response for D = 1 and fs/nQ ^ 10 is different, as 
the force per bead is outside the elastic regime. If we reduce the force per bead for D = 1 
further (filled circles) we obtain finally similar responses as for the larger disks. Note however, 
that linear response requires smaller forces per bead for smaller disks than for larger ones. 

We now consider the mean stresses, i.e. the stress profiles averaged over many realisations 
(different samples and different application points of the force). Far from the source, the self 
averaging discussed above implies that these mean profiles should behave in accordance 
with GET. This is less obvious close to the source, where fluctuations from one realisation 
to the other are large. In Fig. IHl we present the normal mean stresses (axx) and (cTyy) as 
functions of x for different vertical distances y. Similar curves have been obtained for the 
shear stress (cr^y). The agreement with GET (bold lines) is surprisingly good even for small 
distances from the source. It improves further with increasing distance y. Apparently, the 
noise entering in the stress calculation is of (essentially) vanishing mean. While the vertical 
normal stress must have always one peak centered below the source, the horizontal normal 
stress is predicted by classical isotropic theory to show a minimum at x = between two 
peaks for D <^ \y\ <^ h. This is a direct consequence of elasticity. The double peak disappears 
close the fixed surface as there horizontal displacements which cause the tensile horizontal 
forces are suppressed. 

As can be seen from Fig. [7| for the normal stresses, all measured stresses decrease essen- 
tially as the inverse distance from the source (taken aside the expected corrections due the 
finite value of the system size). The two panels given in this figure correspond to measure- 
ments along two straight lines through the source: (a) x = 0, (h) x/y = tan(^^) = ±1. Both 
figures look qualitatively similar. 

More importantly, we compare in both panels both normal average stresses with their 

1 /2 

respective fluctuations from sample to sample Saa/s = {{^ais) ~ i^a/s)^) ■ We note first that 
Saxx ~ Sayy ~ 6axy (the latter relation not being represented in the figure) and that the 
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fluctuations do not depend on the angle 6 of the straight line, but solely on their distance 
r from the source. Surprisingly, we find fluctuations of order of the mean (normal) stresses, 
i.e. the relative fluctuations are of order one close to the source. ^ This striking observation 
is by no means in conflict with the observed self- averaging far from the source, due to the 
different distance dependence of mean stresses and fluctuations. While the former decrease 
(essentially) analytically, our data suggests an exponential decay for the latter. Our fits are 
compatible with a characteristic screening length scale b of order 30. Interestingly, this is 
of same order as the characteristic wave length ^ we have found in j.2jj] for the breakdown 
of the classical eigenmodes. Only for distances somewhat larger than b, the self-averaging 
dominates over the analytical decay of the average stresses, and the relative fluctuations 
vanish eventually. 

In Fig. IHl we discuss in more detail the distribution of stresses along the a; = line through 
the source. Only the vertical normal stresses presented here, as the histograms for 

axx and axy show similar behavior. The normalized histograms have been rescaled and 
plotted versus the natural scaling variable u = cxyy/ {<7yy) which takes out the trivial distance 
dependence of the mean stress. Incidentally, as we know from Fig. IHl we may equally use 
the analytically obtained stress as reference in the scaling variable, without changing the 
reduced histograms. 

Three remarks have to be made here: First, we note that all histograms scale reasonably 
well and the fluctuations, i.e. the width of the unsealed peaks, scale broadly as the mean 
stresses. Closer inspection reveals, however, that the rescaled peak width becomes slightly 
narrower with distance to the source. Both observations are obviously in perfect agreement 
with the previous Fig. [7| where more or less constant relative fluctuation have been found 
due to the large value of the self-averaging length 6 ~ /i. This masks somewhat the different 
functional dependency (analytic versus exponential) of the first two moments of the stress 
distributions. Second, the distributions are more or less symmetric and the mean stress 
corresponds to the maximum of the histogram. This confirms the statement made above 
(Fig. ini), that the fluctuations around the analytical prediction appear to be of vanishing 
mean. Third, although our statistics is certainly insufficient to characterize much better the 
shape of the distributions, specifically the scaling of their tails, a Gaussian distribution can 

^ For symmetry reasons, {<Jxy) = for —> and the corresponding relative fluctuation diverges. 
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be ruled out with the present data. In fact, as shown in the figure, an exponential fit is not 
unreasonable. In this sense the noise is large. 

IV. RESULTS FOR SYSTEMS WITH COMPENSATION FORCES. 

While the previous section was mainly concerned about forces and stresses and their 
distribution, we will now, for the rest of this paper, investigate the displacement of the 
center of mass of the source region. This is the direct route to characterize the local elastic 
properties of an amorphous body. We use here open periodic boundary conditions without 
fixed particles, but with additional small compensation forces on all beads. As before, a 
vertically downwards pointing force acts on source disks of diameter D = 4. 

Fig. ini presents a typical snapshot of the noisy part Su^ = ILs ~ {iLs) °f source dis- 
placements measured in one configuration. For the given box size, the mean displacement 
substracted is roughly four times larger than the average fiuctuation {SuJ^Y^'^ (see Fig. [TT] 
below). Hence, the local elastic properties vary weakly with position. The snapshot (or a 
more detailed histogram) shows the bimodality of the du^ distribution: Few very strong dis- 
placements point downwards in the direction of the force. They are due to some very soft 
spots. The remaining du^ are much smaller and strongly correlated in space. While pointing 
pretty much in all directions, they compensate obviously the net downward component of 
the soft spots. 

We have checked the spatial correlations of the source displacements, by means of the 
(normalized) correlation function (5tt^(r) ■ 5m^(0)) which is summed over all pairs of dis- 
placements of a given configuration, and averaged over the configuration ensemble. In total, 
nearly 4000 responses contribute to the average correlation function presented in Fig. ^1 
For very small distances, the correlation function should become constant since two sources 
of finite disk diameter excite the same beads. Equally, it is expected to become fiat around 
r L/2 due to the periodic boundary conditions. Both limits are in agreement with our 
data. More importantly, the intermediate distance regime D /2 <^ r <^ L/2 may be reason- 
ably fitted in log-log coordinates by a power law slope. Although a somewhat weaker value 
would even fit a larger range of the data we have indicated an exponent —1. Interestingly, the 
same dependency has been observed for the correlation function of the non-affine part of the 
displacement field discussed in Ref . ■ Finally, we remark that an additional non-analytic 
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and possibly exponential regime for 6^r^L/2is conceivable for even larger boxes. This 
is suggested by the existence of the characteristic length b observed in the self-averaging of 
the stresses. Unfortunately, as indicated by the broken line, this is at present not supported 
by our simulations due to the limited system sizes available. 

The reader will have observed that we do not probe here the rigidity of an isolated patch 
of material in a fixed frame. The fluctuations of the source displacements do not depend 
only on the local elements of the dynamical matrix but on a much larger neighborhood, in 
principle the whole system, whose effective size is estimated in the next section. 

V. RESULTS FOR FIXED TOPOLOGICAL LAYERS. 

For solving the linear response directly by means of the dynamical matrix, it is useful 
to renumber and regroup the beads in topological layers around the source disk. All the rii 
beads, interacting directly with the no disk beads, are contained in the first neighbor layer, 
the n2 beads interacting directly with the ni beads of the first neighbor shell are contained 
in the second layer, and so on (there are no direct interactions between the no beads of the 
source and those of the second layer). Both the number of beads of each layer and the mean 
radius R of the fixed particle layer around the source increase linearly with the topological 
rank from the source. The width of a topological layer is of order 3 due to the cut-off of our 
potential and to the weak polydispersity j22|. The last layer of free particles containing ni 
beads, we fix the positions of the — {uq + rii + ... + rii) > remaining beads. 

With this renumbering, the structure of the dynamical matrix becomes more transparent, 
picturing systematically the influence of the subsequent topological layers. While in the last 
section all beads have been allowed to respond to the external load, we study here the 
effect of additional degrees of freedom when more and more topological layers and degrees 
of freedom are allowed to relax. The method used here is a systematic inversion of the 
dynamical matrix. 

As in the last section, we only consider the displacement field at the source. In contrast, 
we use a source containing only one bead {uq = 1). The applied force is arbitrarily set 
to one. Obviously, the response of the source, when all other beads are fixed {I = 0), is 
directly described by the inverse diagonal coefficients of the dynamical matrix, which is a 
function of the local quenched forces and spring constants between neighboring monomers. 
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For / > 0, the displacement U_ of the source can be computed recursively by writing the 
equilibrium equations on all free beads. For 1 = 2, for example, we get after combining all 
the 2 X [rio + ni + 71,2) equilibrium equations in the two directions x and y 

F=(m -M ■ (m -M -M-i -M ) ^ -M ]-U, (1) 

— \ ^=noxno ^=?ioXni y^=nixni ^=mxn2 ^=712 xn2 ^=n2Xni J ^=niX?io / — 

where M.^.^^, is the matrix of size 2ni x 2nj containing all the coefficients of the dynamical 
matrix relating the particles of the layer i with the particles of the layer j. The ratio Us/fs, 
of the vertical component of the source displacement to the applied vertical force, can thus 
be easily computed, with the direct use of the coefficients of the dynamical matrix. 

The vertical component of the source displacement Ug is shown in Fig.^Jas a function of 
the mean diameter 2R of the spherical region around the source which is allowed to respond 
to an external force. The open symbols correspond to the (reduced) mean displacement 
{ug)/ fs, the filled symbols to the fluctuation {6uiy^^ / fg. The squares are for the responses 
measured numerically after relaxation in periodic systems of linear size L = 2R = 104 
without flxed beads fSec. IIV|) . Note that the response at the source depends of course on D 
for all 2R. ^ 

We show, that the mean displacement increases logarithmically with system size (ug) / fs ^ 
0.011 log(2i?/D) in qualitative agreement with continuum theory and this already for sys- 
tems with only one free topological layer around the source (/ = 1). Hence, although / = 
is not sufficient, one can obtain the average local elastic moduli from a surprisingly small 
neighborhood region. Interestingly, the fluctuations level off at much larger distances of or- 
der of 6 = 30 (corresponding to / = 6 topological layers), as can be better seen from the 
inset. (It can be shown that the approach of the large system size limit is exponential.) This 
shows that, at system sizes of the order of the self-averaging length, the fluctuations become 
system size independent. For larger systems, b determines the size of the region responsible 
for the noise in the source displacement fleld. 



^ As different D have been used for the two different boundary conditions, the results from the open periodic 
boundary method have been shifted vertically, taking into account the logarithmic correction suggested 
by continuum theory. 
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VI. CONCLUDING REMARKS. 



We have probed the incremental stress and displacement fields due to a point source 
force acting on two-dimensional amorphous Lennard- Jones solids. Focusing on the linear 
elastic response, this has been done for three different boundary conditions by means of the 
linearized Euler-Lagrange forces (i.e., the dynamical matrix) or by direct minimization of 
the total Hamiltonian. 

We demonstrate that the average stresses and displacement fields compare well with 
the predictions from classical isotropic elasticity, and this already for small distances from 
the source and for small system sizes (Figs. IH| 11111. Contrasting to this, large stress (and, 
hence, strain) fluctuations are found for small distances to the source decreasing (essentially) 
exponentially with distance (Fig. [7j). A surprisingly large length scale 6 ~ 30 is associated 
with this self- averaging with distance. Similarly, the fiuctuations of the source displacement 
fields are found to become system size independent for L = 2R ^ 30. The self-averaging 
length b is of the same order as the characteristic length scale ^ associated with the non-affine 
displacement field under macroscopic strain setting the lower bound for allowing classical 
eigenfrequency calculation to be apphcable. We believe that both length scales express the 
same physical fact and are indeed (up to prefactors) identical quantities. This explains why 
vibration modes corresponding to wavelengths larger than b (and, hence, averaging over 
arger distances) follow continuum theory, while modes with smaller wavelengths do not 



As it has been pointed out, better statistics and much larger box sizes L/2 ^ 6 ~ ,^ would 
be required to establish unambiguously the scaling of the various distribution functions 
discussed here. Especially, an improved correlation function of the source displacements 
would be highly interesting to test the range of the observed power law (Fig. [TUI) . We 
strongly expect the existence of final cut-off at a characteristic length of order b. Additional 
theoretical guidance is also required to explain how such a large length scale arises for 



the fluctuations given the short-range correlation of the dynamical matrix elements 2^|. 
These results should be compared with experiments of nanoscale indentation on glasses 2l| , 
allowing a flrst experimental evidence of large scale fluctuations of the elastic properties in 
amorphous materials. The study of the pressure dependence of our results is currently in 
progress. 
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Coming back to the static properties of granular aggregates (composed of hard, cohe- 
sionless grains in frictional contact) invoked in the introduction, this work suggests several 
computer experiments. As a first step, one should compare the average incremental stress 
and displacement responses. It may be interesting to verify if the double peak structure 
found for the horizontal normal stress is also present in granular matter although there the 
total normal stresses may not become tensile. (A priori this is possible since the incremental 
stresses are perfectly entitled to become negative as long as Coulomb's criterion is not vi- 
olated.) More importantly, the self-averaging properties of granular systems should be put 
to a test. Various experimental facts (especially for forces in vertical columns and silos) Q] 
suggest much larger fluctuations with much weaker self-averaging properties compared to 
amorphous elastic bodies. Finally, it is a matter of debate, if the response to an arbitrary 
weak source does correspond to a Green's function in a strict mathematical sense. Additivity, 
linearity and reversibility of the responses should be tested directly. As the force network 
is subject to incessant restructuring due to the missing permanent grain contacts it may 
not be possible to describe — even in the hydrodynamic limit — the total charging of the 
packing /.near operation Q. 
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APPENDIX A: TWO-DIMENSIONAL ELASTIC RESPONSE TO A POINT 
FORCE. 

At equilibrium the stress state of a two-dimensional elastic material must satisfy the 
(two) force balance equations VjO"ij = F6{x)6{y) with i,j = x,y and F being the external 
point force applied in (0,0). The main assumption of classical elasticity is that the 
system may be entirely described by the continuous displacement field u{x, y). Due to global 
translational and rotational invariance, only the strain field — by definition the symmetric 
part of the gradient of the displacement field — appears in the equations. Moreover in linear 
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elasticity, the stress tensor is related to the strain tensor eij through Hooke's law. Only 
two phenomenological parameters are required for isotropic homogeneous systems, E and u 
(or A and /i). Supposing this, it follows the compatibility equation [2^ d^xC^yy — ^dyyOyy + 
dyyCTxx — J^dxxO'xx = 2(1 + i^)dxy<Jxy and thus, combined with the force balance equations, the 
well-known Laplace equation A^a^x + o-yy) = for {x,y) 7^ (0,0). 

As a specific example we present here the calculation for the point source problem be- 
tween two fixed horizontal walls posed in Fig. ^b). Hence, the displacement field u must 
vanish for \y\ = h. Periodicity and symmetry of the simulation box in horizontal direction 
impose a solution periodic and symmetric (odd or even) in x. These boundary and symmetry 
conditions can be readily reformulated in terms of the stresses. 

To obtain the elastic response, the idea is to use the method presented in ^ in the case of 
an elastic layer submitted to a force localized at its surface. We divide our medium into two 
parts: part 1 above, part 2 below the point source. The continuity of the displacement field 
along the fictitious dividing line requires u^x\x,0) = u^x\x,0) and u^y\x,0) = u^\x,0). 
The point force is taken into account by imposing (Jy^Jix, 0) = cry^{x, 0) — p{x) where p{x) 
is the vertical external pressure. An additional constraint is imposed by the continuity of 
the shear stress a2y{x, 0) = cr^ij{x, 0). Note that the continuity of Ux ai y = together with 
the discontinuity of cTyy imposes the discontinuity of axx- Imposing a continuous axx would 
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yield to a discontinuous displacement field Ux with large scale vortices 

The stress tensor components are decomposed into a base of harmonic functions, typically 
affine functions or product of trigonometric functions with exponentials. Taking into account 
the X < — > —X symmetry, we look for a solution of the type 

+00 

^i'/^ + 4''^ = $^cos(gx) (a(^'2)e^^ + b'^'''^e"^y) 



n=0 

+00 



- ^S''^ = $^cos(gx) {qy {a^^'^^e'^y - b^^'^^e-'^y) + 2 {c^^'^^e'^y - d^''^^ e^^y)) 

n=0 
+00 

= ^sin(gx) {qy/2 [a^^'^'^e'^y + b^^'^^e-'^y) + {c^^'^'^e'^y + d^^'^^e-^y)) . (Al) 

n=0 

^(1,2)^^(1,2)^^(1,2) g^^^ ^(1,2) 

are coefficient functions depending on the frequency q = nAq. The 
latter is quantized with Aq = 27r/L, due to the finite horizontal width L of the layer, and 
periodic boundary conditions. A similar looking ansatz for the corresponding displacement 
fields is readily obtained. 
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The vertical external pressure is expressed in the same form, as 



+00 



p{x) = cos(gx)s(g).Ag 



n=0 



with p{x) either equal to a gaussian ^s(g) = F/Tre*^""^^"^/^-*^, or a uniform force of width 2a 
{s{q) = F sm{qa) / (naq)) . 

Replacing the general expressions for the stresses and the displacement fields in the 
equations characterizing the boundary conditions and the continuity at the fictitious dividing 
line, leads to the following explicit expressions for the eight functions a*-^'^-*, b^^''^\ c^^'^-* and 
(^(1:2) (depending on s{q), qh and u 

s{q)Aq ((1 + u).2qh + 3 - u) e'^'^^ + (z/ - 3)e-^i'' 



a\q) 
a\q) 

c\q) 
c\q) 



4 D{q) 
s{q)Aq ((1 + iy).2qh - 3 + z/) e-^^^' - (^^ - 3) 

s{q)Aq (-(1 + iy).2qh + 3 - e'^"'^ + (^^ - 3) 
~^ 

s{q)Aq ((1 + u).2qh + 3 - e^^^^ + (z/ - 3)e-^«'* 
~^ W) 

s{q)Aq (2(1 + u)\^h^ + (1 - u^)2qh + (1 - z/)(3 - u)) e-^^^ + (1 - - 3)6"^^^ 
8(1 + 1%) 
rf2(g) 

g(g)Ag (2(1 + ufq^h^ - (1 - z/^)2g/i + (1 - z/)(3 - u)) e"^^^ + (1 - z/)(z/ - 3) 
8(1 + z/) D{q) 
= d\q) 

w.thZ.(,).2,/.e--^ + ^e-^^ 

Substituting these coefficient functions back into the general ansatz for stress and dis- 
placement fields, one obtains explicit expressions for the stress and the displacement fields. 
We have drawn numerically these expressions to get the theoretical fits presented in Sec. IIIII 
Note that unlike the boundary condition studied in \}\ , the coefficient functions depend now 
on u. The displacement field u is proportional to 1/E. The solution also depends on the 
system height 2h, on the size a of the source, and on the width L, the latter due to the 
quantization of the Fourier integration. 
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FIG. 1: (Color online) Sketch of two boundary conditions of interest for measuring the response 
to an additional point force source (a) The source may be applied to the free upper surface of a 
prestressed aggregate formed at constant gravity on a rigid bottom plate (possibly containing some 
stress transducers). This setup has been studied extensively recently {^,0,01 in order to determine 
the static response of packings of (hard and cohesionless) granular matter, (b) One of the three 
boundary conditions studied in this paper. The source is applied within a macroscopically isotropic 
and homogeneous "computer solid" in a periodic simulation box of linear size L. The center of the 
source defines the origin of the coordinate system {x,y). For mechanical stability we either apply 
a compensation force of —f^/N (N being the total number of beads) to all particles or we freeze 
some particles (gray beads) as shown in the right panel. We study the response of amorphous 
packings of carefully quenched (slightly polydisperse) Lennard- Jones beads. Obviously, this is a 
further simplification with regard to the granular material case with its more intricate non-linear 
(static friction) particle interactions. 
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FIG. 2: (Color online) Snapshot of the incremental forces in one small periodic box containing 
N = 1000 particles generated by the source applied on all the beads within the disk indicated. We 
have chosen here a disk diameter of 4 particle sizes. The line width between interacting beads is 
proportional to the incremental forces. (Only forces larger than 0.02 have been drawn for clarity.) 
Black (gray) lines correspond to incremental compressive (tensile) stresses. Also indicated on top 
and bottom are the beads fixed to balance the source. The snapshot shows that the forces generated 
by one additional source are strongly heterogeneous and resemble qualitatively the "force chains" 
known from granular matter ^. 
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FIG. 3: (Color online) Snapshot of the (reduced) displacement field 6u = u — Ucet generated in 
the same configuration as in Fig. |2 We have substracted from the total displacement field u, the 
displacement field Ucet obtained analytically from classical continuum elasticity theory (GET). The 
difference from continuum theory is quite marked for the displacement of beads on the "force 
chains" of Fig. |2 On larger distances rotatory structures become visible — quite similar to the 
ones obtained from the non-affine part of the displacement fields under macroscopic strain 2^ . 
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FIG. 4: (Color online) LEFT: Snapshot of the (reduced) displacement field 5u = u — Ucet super- 
imposed with the noise component of the incremental stresses (that is cr — (J quenched ~ (^cet)- The 
chosen configuration is the same as in Fig. |2l and Fig. |21 Stresses are represented by ellipses whose 
large principal axis is proportional to the largest eigenvalue of the local (incremental) stress tensor. 
The small axis is proportional to the smallest eigenvalue of this stress tensor. The directions of the 
axes of the ellipses give the main directions of stress. The arrows represent the displacement field, 
as in Fig. 01 

RIGHT: Histogram of the angles 9 between the local (reduced) displacement, and the main direction 
of the (incremental) stress tensor. The histogram is peaked around zero, with a broad distribution 
oc 6. It has been obtained from 10 configurations of = 10000 particles. 
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FIG. 5: (Color online) Unaveraged vertical stress Uyy close to the bottom plate {y = 50) caused 
by sources of various disk diameters D (as indicated in the figure) at the same position of one 
configuration of linear size L = 104. The total vertical stress has been used as a normalization. The 
open symbols correspond to a total applied force fs = 10, the filled circles are for a source with 
D = riQ = 1 and fs/riQ = 1. The bold line shows the theoretical prediction. It corresponds also to 
the statistical average (see Fig. E|. The linear responses for D = 4,6 and 8 are perfectly identical. 
This only applies as long as the force per bead remains sufficiently small. 
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FIG. 6: (Color online) Averaged incremental stresses axx (left) and ayy (right) versus x for different 
vertical distances y from the source. The boundary conditions indicated in Fig. 1(b) are used. 
Data from configurations containing N = 10, 000 beads in boxes of L = 104 is averaged over 220 
independent measurements and compared with the predictions from classical elasticity (bold lines). 
The agreement is surprisingly good even for small y and improves systematically with increasing 
source distance. 
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FIG. 7: (Color online) Comparison of the (incremental) stress fluctuations daaf3 = 
/ \ 1/2 

ii'^afs) ~ i'^afs)'^] with the mean vertical normal stress (cyy). This is done in two directions 
through the source: (a) along the vertical line {x = 0) and (b) for \x/y\ = 1. We note that mean 
stresses and their fluctuations scale quite differently with distance r from the source in both direc- 
tions. For small distances we find relative fluctuations (^Cq,^/(cJq,^) of order one. While the mean 
stresses decrease, as expected in 2D, essentially as 1/r, the fluctuations are found to be well fitted 
by an exponential decay exp(— r/6) with b w 30. 
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FIG. 8: (Color online) Normalized distributions of incremental vertical normal stresses for different 
distances y > 0, along the vertical line through the source {x = 0). The histograms are plotted 
versus u = cryy/{ayy). The data is averaged over 220 independent point source experiments. The 
histograms are more or less symmetric. With increasing distance, the rescaled distributions become 
systematically narrower, in agreement with the previous figure, however, the effect is weak. For 
y < 40, the tails of the histograms are not Gaussian, but roughly exponential as indicated by the 
bold lines. 



26 




FIG. 9: (Color online) Snapshot of the reduced displacement vectors 6ug = Ug — (Ug) of the center 
of mass of the source region. (The size of the arrows is proportional to the length of the reduced 
displacement vector.) This data have been obtained for completely periodic boundary conditions 
where no particles have been fixed. The vector field varies greatly in size and direction. Closer 
inspection shows strong spatial correlations. 
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FIG. 10: (Color online) Spatial correlation function {Sus{r)-5ug{0)) / {Ug"^) of the source displacement 
vector SUg = Ug — {Us) with r being the distance between source terms within the same configuration. 
Note that the correlation function is normalized. The average is taken over a total number of 
nearly 4000 linear responses using open periodic boundary conditions without fixed particles. As 
indicated by the bold line, the correlation function decreases essentially like the inverse distance 
for D/2 <^ r <^ L/2. It becomes constant for smaller and larger distances. We strongly suspect an 
additional exponential cut-off (dashed line), however, our data is too noisy and, more importantly, 
L is too small to demonstrate this unambiguously. 
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FIG. 11: (Color online) Vertical component of the displacement of the source center of mass Ug 
for a given total force fg as a function of system size 2R. The open symbols correspond to the 
mean displacement (n^)^/^//^, the filled symbols to the fluctuation ((5n^)^/^//s. The circles are 
for the boundary condition with fixed topological layers discussed in Sec. El R being the mean 
distance to the fixed border shell. The squares are for responses in periodic systems of linear size 
L = 2R = 104 without fixed beads (Sec. IIV|) . The mean displacement increase logarithmically 
with system size {ug)/ fs ~ 0.011 log(2i2/L') in agreement with theory. The fluctuations level off at 
distances of order oih = 30, as can be better seen from the inset. 
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